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LONG-TERM  GOALS 

For  an  V-element  array  (Fig.  1(a)),  methods  such  as  beamforming  and  singular  value  decomposition 
rely  on  estimation  of  the  sample  covariance  matrix,  computed  from  M  independent  data  snapshots.  As 
M  — >  oo ,  the  sample  covariance  is  a  consistent  estimator  of  the  true  population  covariance.  However, 
this  ideal  condition  cannot  be  met  in  most  practical  situations,  '  in  which  large-aperture  arrays  operate 
in  the  presence  of  fast  maneuvering  interferers,  or  with  towed/drifting  arrays  strongly  affected  by 
defonnation  or  array-depth  perturbations.  The  long-term  goal  of  this  effort  is  the  development  of 
physically  motivated  models  to  statistically  describe  the  eigenstructure  (eigenvalues  and  eigenvectors) 

of  sample  covariance  matrices  in  sample-starved  settings,  and  the  use  of  those  models  for  performance 

analysis  and  improvement  of  array  processing  methods.  To  this  end,  mathematical  tools  developed  in 
the  context  of  Random  Matrix  Theory  (RMT)3"  (mostly  focused  in  the  regime  N~M)  and  High 

IQ. 

Dimension,  Low  Sample  Size  (HDLSS)  array  processing  "  (which  considers  N»M)  are  applied  to 
obtain  statistical  descriptions  of  sample  eigenvalues/eigenvectors  and  how  those  quantities  differ  from 
the  (true)  population  eigenpairs.  Additional  long-term  goals  are  exploiting  the  information  carried  by 
sample  eigenvectors  for  the  improvement  of  estimators  of  the  sample  covariance  matrix  (i.e.,  signal 
versus  noise  subspaces),  and  for  quantifying  local  stationary  in  array  data  (Fig.l  (b)). 

OBJECTIVES 

•  Developing  beamforming  techniques  for  snapshot-deficient  scenarios  with  moving  targets  such  as 
the  one  illustrated  in  Fig.  1 .  Unlike  previous  research,  this  work  incorporates  asymptotic  results  " 
for  the  sample  eigenvectors  and  how  they  deviate  from  their  true  population  eigenvectors. 

•  Developing  multi-hypothesis  tests  for  estimation  of  parameters  such  as  azimuth  and  power  of 
targets  in  the  watercolumn,  number  of  targets,  and  their  speed,  based  on  asymptotic  limits  for  the 
sample  eigenvalues,  sample  eigenvectors,  and  eigenspace  projections  for  N»M.  A  key 
component  of  this  research  is  the  analysis  of  the  impact  of  source  movement  (i.e.,  time-dependent 
position)  on  the  eigenspaces  associated  to  signal  and  background  noise. 
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Figure  1  (a)  Experimental  scenario  considered  in  this  research,9  consisting  of  multiple  maneuvering 
sources  detected  by  an  N-element  array  after  collecting  M  snapshots.  This  work  considers  array 
processing  tools  for  cases  in  which  N»M.  (b)  Data  snapshots  are  collected  at  discrete  times  ta,  ta+1, 
etc.  Data  local  stationarity  occurs  at  time  intervals  of  little  or  no  variation  of  target  azimuth  (e.g., 

ta<t<ta+2  and  th<t<th+s). 


•  Developing  statistical  methods  for  data  segmentation  into  intervals  with  local  stationarity9'10  as 
shown  in  Fig.  1(b),  capable  of  distinguishing  between  true  variations  in  the  covariance  structure 
(due  to  underlying  time-varying  statistics),  as  opposed  to  variations  introduced  by  the  lack  of  data 
snapshots.  Since  estimation  of  the  sample  covariance  matrix  assumes  local  stationarity  of  the  data, 
this  segmentation  approach  will  provide  a  bound  for  the  largest  number  of  snapshots  to  be 
included  while  estimating  covariances. 

•  Obtaining  an  eigenvector-based  estimator  of  the  signal  subspace11  to  improve  the  performance  of 
subspace  array  processing  methods.  In  previous  ONR-sponsored  work,  determining  the  signal 
subspace  has  been  carried  out  by  defining  eigenvalue-based  rank  estimators.1’  In  contrast,  the 
work  presented  here  considers  the  benefits  of  including  the  infonnation  carried  by  sample 
eigenvectors. 


APPROACH 

1)  Defining  a  snapshot  model:1  Let 
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be  a  data  snapshot  at  time  tm  at  an  A-elcmcnt  array,  where  Wm  is  the  Nxl  vector  corresponding  to  the 
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background  noise  of  power  level  crw  ,  yqm  is  a  complex  zero-mean  Gaussian  distributed  random 
variable  with  variance  crq  representing  the  source  level  of  the  qh  target  in  the  watercolumn  (see 
Fig.  1(a)),  and  V qm  is  the  azimuth-dependent  replica  vector  for  the  q'h  target,  defined  as 
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where  CO  is  the  frequency  in  radians/s,  d  is  the  array  inter-element  spacing,  and  co  is  the  water  sound 
speed.  The  population  covariance  at  time  ta  is 
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Using  eq.(3),  data  realizations  (i.e.,  M  snapshots)  can  be  generated  from  complex  normally  distributed 
vectors  z  as  xm  =  R'/2  z  .  The  simulated  data  obtained  by  this  approach  is  used  to  study  the  behavior  of 
the  sample  covariance  matrix,  estimated  as: 
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From  eq.  (3),  the  relative  impact  of  Q ,  Q  ,  a2  /  a2 ,  N  ,  and  M  on  the  behavior  of  sample  eigenvalues 

( X  j. . .  A  N  )  and  sample  eigenvectors  ( u  v . .  u  N  )  will  be  analyzed  by  Monte  Carlo  studies  and  contrasted 
to  predictions  from  RMT3'6  and  HDLSS. 7-8 

2)  Quantifying  eigenspace  dynamics:9'10’11  eigenspaces  are  defined  as  the  vector  space  spanned  by 
Q  eigenvectors.  Signal  stationarity  can  be  quantified  by  tracking  time-dependent  variations  in  the 
eigenspace  defined  by  target-related  population  eigenvectors,  as  illustrated  in  Fig.  2.  In  this  example 
the  decision  of  data  stationarity  is  based  on  whether  the  population  eigenvectors  have  changed 
direction  over  time  (i.e.,  u  la=  u  lb  and  u  2a=  u  2b )  or  not.  Since  population  eigenvectors  are  not 
observable  variables,  an  eigenspace  distance  based  on  sample  eigenvectors  is  defined  as:9'10 
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which  has  bounds  0  <  rab  <  1  since  the  trace  operator  is  bounded  between  0  and  1 .  rab  =  0  occurs  when 
the  eigenspaces  at  times  ta  and  th  are  perfectly  lined  up  indicating  data  stationarity.  On  the  other  hand, 
rab  =  1  results  when  such  eigenspaces  are  orthogonal  to  each  other,  indicating  data  non-stationarity.  A 
statistical  decision  rule  for  data  stationarity  can  be  stated  as  the  complementary  hypothesis 

H0 :  stationary  data  — »  -L  <  rab  <  +L 

>  (6) 

Hl  :  non  stationary  data  ~^-L>  rab  or  rab  >  +L 

where  the  stationarity  bounds  ±L  =  jur  ±  <jr  are  a  function  of  the  theoretical  mean  ( /ur )  and  standard 

on 

deviation  ( cr,. )  of  rah .  Based  on  theoretical  results  from  HDLSS  '  it  is  possible  to  obtain  theoretical 
expressions  for  these  quantities  on  the  assumption  of  data  stationarity  (hypothesis  H0 ). 
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Figure  2  Illustration  of  the  data  stationarity  test  proposed  in  this  work,9'10  which  consists  on 
tracking  variations  in  the  direction  of  population  eigenvectors  at  times  ta  and  th.  Since  population 
eigenvectors  are  not  observable  quantities,  eigenspace  variability  can  be  tracked  using  the 
eigenspace  distance  metric  in  eq.(5),  based  on  observable  sample  eigenvectors. 


3)  Obtaining  an  eigenvector-based  estimator  for  signal  and  noise  subspaces:10  Subspace 
beamforming  techniques  rely  on  projecting  noisy  data  snapshots  into  reduced  rank  signal  eigenspaces. 
In  this  work,  we  consider  the  statistic  properties  of  entirely  random  eigenspaces  to  identify  specific 
eigenvectors  with  signal-bearing  information.  To  explain  this  concept,  consider  the  following  Monte 
Carlo  experiment  illustrated  in  Fig.  3:  given  an  /V-dimcnsional  sphere,  randomly  select  K  vectors  and 
measure  the  angle  a  between  these  vectors  and  an  arbitrary  reference  vector  (shown  in  red).  It  can  be 
shown10’12'14  that  for  random  vectors  (i.e.,  vectors  with  no  information  content),  the  angle 
«min  =min(«j ,...,aK)  has  a  probability  distribution  function10 
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This  concept  can  be  extended  to  K  eigenvectors  estimated  from  data  snapshots  collected  at  an  /V-scnsor 
array:  sample  eigenvectors  corresponding  to  noise-only  data  exhibit  the  same  statistics  as  eq.(7),  while 
signal-bearing  eigenvectors  result  into  left-skewed  distributions,  as  illustrated  in  Fig.  3.  Therefore, 
eigenvectors  with  amin  lower  than  a  user-defined  threshold  should  be  considered  as  potentially 
“informative”.  In  this  work,  we  propose  an  algorithm  described  below  that  utilized  this  criteria  for 
estimation  of  the  signal  subspace  prior  to  beamforming. 
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Figure  3  Example  of  the  approach  proposed  in  this  work  to  identify  signal-bearing  vectors:10  the 
minimum  angle  between  a  set  of  K=3  noise-only  N-dimensional  vectors  and  a  reference  vector  (red) 
exhibits  different  behavior  than  the  angle  resulting  from  signal-bearing  vectors  (i.e.,  those  vectors 
with  structure  that  resemble  the  reference).  In  this  research,  this  concept  is  extended  to  sample 
eigenvectors  as  a  tool  to  estimate  the  signal  subspace,  as  well  as  its  complementary  noise  subspace. 

WORK  COMPLETED 

1)  An  eigenvector-based  metric  to  quantify  data  stationarity9  (based  on  the  concept  illustrated  in 
Fig.2)  was  proposed  for  an  experimental  scenario  consisting  of  a  single  dominant  interferer  and 

multiple  quiet  targets  of  interest.  A  generalization  of  this  approach  to  the  case  of  Q  dominant 
interferers  has  also  been  proposed10  and  the  theoretical  results  have  been  applied  to  improving  a 
reduced  rank  beamformer  and  applied  to  simulated  and  experimental  data. 

2)  A  sample  eigenvector-based  estimator  for  the  signal  subspace  has  been  proposed10  based  on  the 
idea  illustrated  in  Fig.  3.  Beamforming  results  with  simulated  data  illustrate  application  of  this 
theoretical  result  to  improve  target  detection  in  snapshot-starved  scenarios. 

Examples  of  results  corresponding  to  applications  #1  and  #2  are  shown  in  the  next  section. 

RESULTS 

The  array  processing  techniques  developed  in  this  research  have  application  to  the  improvement  of 
subspace  beamfonning  processors,  in  which  the  data  snapshots  are  projected  into  particular  subspaces 
of  interest  such  as 

y m  =  p<  )  X,„  »  (8) 

where  P(  ’  is  a  subspace  projection  operator  obtained  from  sample  or  population  eigenvectors,  defined 
in  Table  1  for  several  cases  of  study. 
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Table  1:  Definition  of  three  data  projectors  P'  and  adaptive  beamforming  weight  vectors  w°(#) 
used  in  this  report  for  comparison  of  array  processing  techniques. 


Projector 

Adaptive  weights 

Description 

n= 1 
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1  a+M-\  r 
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Projection  into  true  signal 
eigenspace.  Used  as 
benchmark  for  simulations 

Pw  =  f>  „(Mo)uf(Mo) 
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i  a+M0- 1  r  .. 

R  _  1  V  4-inMGC/1° 
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Projection  into  rank  Q 
subspace.  Sample 
eigenvectors  estimated 
from  Mo  snapshots. 
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^  m=a 

Projection  into  rank  Q 
subspace.  Orthonormal 
basis  u„  estimated  by  the 
algorithm  described  below. 

Once  the  adaptive  weights  have  been  obtained  (see  Table  1),  beamforming  results  have  the  form 

.  (?) 

This  report  discusses  two  beamforming  approaches  that  exploit  either  data  stationarity  or  signal 
information  content  in  sample  eigenvectors: 

1)  Stationarity-aided  beamforming:9'10  Application  of  the  stationarity  metric  in  eq.(5)  for 
improvement  on  the  estimation  of  adaptive  beamfonning  weights  is  demonstrated  in  this  section  with 
simulated  and  experimental  data.  The  simulated  setting  allows  to  compare  the  benchmark  beamformer 
BFffe)  to  BF(MMa°=M) ,  which  is  the  traditional  MVDR  adaptive  beamformer  with  white  noise  gain 

constraint.  This  is  also  contrasted  to  BF^°~M)  in  which  M0  is  detennined  according  to  the  data 
stationarity  criterion  in  eq.(6).  The  simulated  example  consists  of  two  loud  interferers  (20  dB  and  25 
dB)  and  four  quiet  sources  with  powers  ranging  between  2  dB  and  5  dB  above  the  noise  level  of  0  dB. 
Figure  3  (a)  shows  target  azimuths  as  a  function  of  time.  The  data  stationarity  criterion  in  Fig.3(b)  is 

dominated  by  the  loud  interferers,  so  Q-  2  in  eq.(5).  This  is  evident  for  0<t<l  30  s,  for  which  both 
interferers  have  constant  azimuth,  resulting  in  a  large  interval  of  data  stationarity.  Another  interesting 
feature  is  observed  around  450  s,  where  both  interferers  are  undistinguishable  as  a  result  of  being 
located  at  the  same  azimuth.  In  this  case,  Q  =  1  and  eq.(5)  predicts  accurately  a  sharp  drop  of  the 
stationarity  metric.  Figure  3(c)  shows  the  benchmark  beam  former  BFf‘‘e> ,  in  which  all  targets  can  be 
clearly  visualized.  Figure  3(d)  shows  BFffla°=M) ,  obtained  by  estimating  the  adaptive  weights  using 

M=10  snapshots  at  each  time.  Compared  to  panel  (c),  noise  levels  are  higher  and  target  angular 
resolution  is  lower.  Figure  3(e)  shows  the  improved  results  obtained  by  estimating  beamforming 
weights  using  Ma  >  10  snapshots  according  to  the  data  stationarity  bounds  in  Fig.3  (b).  The  results 
resemble  those  in  panel  (c),  with  similar  background  noise  level  and  angular  resolution. 
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Figure  3  Example  of  eigenspace  stationarity  in  simulated  data:10  (a)  Time-dependent  azimuth  of  six 
maneuvering  targets  (red=loud  sources);  (b)  Stationarity  metric  rab  (solid  black  line)  compared  to 
theoretical  bounds  (dark  grey  box).  Horizontal  boxes  indicates  segments  of  data  deemed  as 
stationary  according  to  the  criterion  in  eq.(6);  (c)  BFffe]  (d)  BFff=M  1  with  Ma  =10  snapshots; 

(e)  BF^"-m ] ,  with  M0  >  1 0  snapshots  according  to  the  stationarity  segments  in  (b). 


The  stationarity-aided  beamforming  approach  was  also  applied  to  experimental  data  from  the  Cal  Ops 
portion  of  the  Shallow  Water  Array  Performance  (SWAP)  experiment.14  Although  a  towed  active 
source  was  utilized  during  CalOps  experiment,14  analysis  of  the  full  data  set  is  still  ongoing  and  in  this 
report  we  only  consider  data  segments  with  sources  of  opportunity  (shipping  traffic).  Beamforming 
results  for  data  collected  on  September  7th,  2007  are  shown  in  Fig.  4.  The  data  was  collected  on  a  100- 
sensor  horizontal  line  sub-array,  for  a  total  175  m  aperture.  Figure  4(a)  shows  the  application  of  the 
stationarity  metric  to  the  dominant  time -varying  eigenspace.  Similar  to  the  previous  simulated 
example,  rab  exhibits  variability  due  to  the  lack  of  snapshots  used  to  compute  the  eigenspaces  at  each 
time,  as  well  as  due  to  true  variations  in  the  underlying  data  statistics  from  surface  ship  movement. 
Figure  4  (b)  shows  BF^°=M)  for  M0  =  5  snapshots,  allowing  the  detection  of  multiple  targets  in  the 

watercolumn.  Improved  results  BF^°~M)  are  shown  Fig.4(c),  which  exhibit  reduced  background  noise 
levels  and  better  detection  capabilities  of  quiet  targets  than  panel  (b). 


7 


-20  -10 


?  dB 


Figure  4  Eigenspace  stationarity10  applied  to  experimental  data  from  the  SWAP  experiment.14  The 
data  was  collected  at  a  100-sensor  horizontal  line  sub-array,  at  415  Hz.  (a)  Stationarity  metric  rab 
(solid  black  line)  compared  to  theoretical  bounds  (dark  grey  box).  Horizontal  boxes  indicates 
segments  of  data  deemed  as  stationary;  (b)  BFf4"=M '  with  Ma  =  5  snapshots;  (c)  BF^°~M) ,  with 
Ma>  5  snapshots  according  to  the  stationarity  segments  in  (a). 


2)  Eigenvector-based  beamforming  ( EVB): 10  as  illustrated  in  Fig.3,  eigenvectors  with  information 
content  related  to  targets  in  the  watercolumn  exhibit  different  statistical  properties  to  those  related  to 
background  noise.  In  this  research  we  propose  an  eight-step  algorithm  that  exploits  this  infonnation 
content  to  estimate  the  signal  subspace  for  data  projection  prior  to  beamfonning: 


Estimating  the  eigenbasis  uj  required  for  the  proposed  EVB  beamformer: 

Step  1:  Given M data  snapshots,  compute  the  sample  eigenvectors  u u  w  . 

Step  2:  Compute  the  steering  vectors  \(0  f )  for  a  fine  grid  of  steering  angles  -90°  <0 f<  90° . 

TF 

Step  3 :  Define  TF  as  the  solution  of  F  =  j  fa  (ocmm)damm  ,  where  Fisa  user-defined  constant 

0  <  F  <  1 .  Notice  that  F  is  ultimately  related  to  the  percentage  of  false  detections  at  the  output  of  the 
proposed  beamformer.  For  example,  F=0  results  in  0  informative  eigenvectors  (i.e.,  no  false  alarms, 
but  also  misses  all  true  detections),  while  F=1  indicates  that  all  eigenvectors  are  informative,  leading 
to  large  number  of  false  alanns. 

Step  4:  Initialize  sk  =  u*.  for  k  =  1..M  . 

Step  5:  Find  0f  for  which  acos(|  sk\(0f)  \  )<TF. 
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Step  6:  Remove  \(Of)  from  s^as:  s*.  =  si-|sfv(#/)|  \{6f)  for  k  =  \..M  . 

Step  7 :  Repeat  from  step  5  until  acos(|  sf  \{6 1 )  |  )>  TF  for  all  0 f . 

Step  8:  The  resulting  non-orthogonal  vectors  gk  =  u/;  -s^  (where  k  =  1..M),  span  the  signal 
eigenspace.  Applying  singular  value  decomposition  to  g*  gives  the  orthonormal  basis  ,  where 
k  =  I...Q  and  Q  is  automatically  detennined  by  the  column  rank  of  g# . 


Notice  that  the  main  idea  behind  this  algorithm  is  the  removal  of  plane -wave  wavefronts  from  the 
sample  eigenvectors  (step  6).  Selecting  which  wavefronts  are  to  be  removed  is  done  based  on  the 
criterion  in  step  5.  By  step  8,  st  spans  the  noise  subspace  while  gk  spans  the  signal  subspace. 

The  eigenvector-based  beamformer  described  in  Table  1  was  applied  to  simulated  data  by  projecting 
the  data  into  the  uj  basis  estimated  by  the  previous  algorithm.  An  example  of  this  result  is  shown  in 

Fig.5,  in  which  120  Monte  Carlo  beamfonning  realizations  were  obtained  for  an  experimental  scenario 
consisting  of  a  100-sensor  horizontal  line  array  in  the  presence  of  6  sources  with  azimuth 
0°,1.3o,3o,5.5o,7°,  3°  ,  and  8° .  Each  Monte  Carlo  realization  consists  of  generating  M=10  data 
snapshots  and  applying  adaptive  beamfonning.  Figure  5  (a)  shows  BF^°=M)  (i.e.,  the  traditional 
MVDR  result).  Figure  5  (b)  shows  the  proposed  BF{U faVB)  obtained  by  using  a  threshold  F  =  0.05  (see 
step  3  in  previous  algorithm).  A  detailed  view  of  the  beamfonning  results  is  shown  in  Fig. 5(c)  for  a 
single  Monte  Carlo  realization.  In  this  case,  BF^aVB)  (red  line)  accurately  detects  five  out  of  six  targets, 
missing  the  one  at  1 .3° .  In  addition,  an  evident  false  detection  can  be  seen  at  -  43° .  However,  BF^faVB) 
yields  sharper  detections  that  allow  clear  view  of  the  number  of  sources  used  in  the  simulation. 

Ongoing  research  on  the  proposed  eigenvector-based  beamfonner  is  on  relating  the  user-selected  value 
of  F  to  the  percentage  of  false  detections  at  the  beamfonner  output,  in  such  a  way  to  provide  a 
guaranteed  false  alarm  rate. 

Figure  6  shows10  a  comparison  between  MVDR(Fig.6  (a))  and  EVB  (Fig.6  (b)  and  (c))  beamformers 
applied  to  experimental  data  from  the  SWAP  experiment.14  The  results  for  EVB  were  computed  using 
two  false  detection  levels:  with  F  =  0.05  in  (b)  only  the  most  prominent  targets  are  visualized,  yielding 
fewer  false  detections  compared  to  BF^°=M)  in  (a).  As  in  the  simulated  case,  target  detections  are 

sharp,  suggesting  an  improved  azimuth  resolution.  By  increasing  the  false  detection  level  to  F  =  0.1 , 
the  number  of  target-related  detections  increases  at  the  expense  of  a  higher  number  of  false  peaks. 
However,  even  in  this  case  EVB  significantly  reduces  the  number  of  false  detections  when  compared 
to  MVDR  in  Fig.6  (a).  Details  of  the  improvement  in  azimuth  resolution  of  the  EVB  beamfonning  are 
shown  Fig. 7,  which  conesponds  to  a  zoom-in  view  of  the  results  in  Fig.6. 

IMPACT/APPLICATIONS 

This  work  will  benefit  adaptive  beamfonning  techniques,  in  particular  algorithms  that  rely  on  accurate 
estimation  of  covariance  matrices,  matrix  rank,  and  subspaces.  This  research  aims  to  enhancing  the 
sensitivity  to  detect  genuine  variations  (trends)  in  signal  statistics,  by  applying  results  from  a  growing 
body  of  research  on  the  asymptotic  behavior  of  eigenvectors  of  sample  covariance  matrices.3'8  This 
research  provides:  data-driven  bounds  for  data  segmentation  based  on  data  stationarity  and  improved 
estimators  of  signal  vs  noise  subspaces  of  a  sample  covariance  matrix  computed  from  few  snapshots. 
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Figure  5  Improved  beamforming  by  exploiting  sample  eigenvector  signal  content:10  (a)  MVDR 
beamformer  BF^°=M)  using  M  =  10  snapshots/Monte  Carlo  realization;  (b)  Proposed  eigenvector- 

based  beamformer  BF(Jam  ;  (c)  BF^°=M)  (black  line)  versus  BFff11  1  (red  line)  for  the  first  Monte 

Carlo  realization  in  (a)  and  (b).  Horizontal  dashed  lines  indicate  the  azimuth  of  the  six  sources  used 

in  this  simulation. 


Figure  6  Eigenvector-based  (EVB)  beamforming  applied  to  experimental  data  (f=205  Hz)  from  the 
SWAP  experiment  collected  at  a  sub-array  with  N=60  hydrophones:14  (a)  MVDR  beamformer 
BF^°=M)  using  M  =  10  snapshots;  (b)  BFfa>li]  with  false  detection  level  F  =  0.05;  (c)  BFjfam  with 

F  =  0. 1 .  For  comparison  between  beamformers,  results  are  normalized  with  maximum 

corresponding  to  0  dB. 
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Figure  7  Zoom-in  view  of  Fig.  6,  showing  the  improved  azimuth  resolution  obtained  by  the  EVB 
beamformer  in  (b)  and  (c)  compared  to  the  MVDR  beamformer  in  (a).  Over  plotted  lines  indicate 
actual  ship  lanes  obtained  from  AIS  data  collected  during  the  SWAP  experiment. 14 


RELATED  PROJECTS 

None. 
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